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Employing a recently developed method that is numerically accurate within a model space sim- 
ulating the real-time dynamics of few-body systems interacting with macroscopic environmental 
quantum fields, we analyze the full dynamics of an atomic system coupled to a continuum light-field 
with a gapped spectral density. This is a situation encountered, for example, in the radiation field in 
a photonic crystal, whose analysis has been so far been confined to limiting cases due to the lack of 
suitable numerical techniques. We show that both atomic population and coherence dynamics can 
drastically deviate from the results predicted when using the rotating wave approximation, partic- 
ularly in the strong coupling regime. Experimental conditions required to observe these corrections 
are also discussed. 
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I. INTRODUCTION 

In the past decades, there has been an enormous in- 
terest in developing methods to control and modify the 
light field. In most cases, this has been achieved by 
building specific materials where the light matter inter- 
action strongly modifies the characteristics of the free 
electromagnetic field. These new materials result in in- 
teresting technological applications as well as theoreti- 
cal challenges. Important examples are microcavities 
and metamaterials , which are designed respectively to 
strengthen the light-matter interaction, and to tailor the 
refraction index altering the flow of light in a non-trivial 
way ii. 

Within this scenario, photonic crystals (PC) provide 
one of the most interesting examples of artificially engi- 
neered materials. PCs [^-tZj (see @, Q for basic reviews) 
are periodically microstructured compounds which tailor 
the vacuum electromagnetic density of states, creating 
specific frequency ranges, or gaps, where it vanishes. Re- 
duction or suppression of the density of states within the 
band gap facilitates light localization and trapping in a 
bulk material as well as the inhibition of sponta- 

neous emission over a broad frequency range. Addition- 
ally, the density of states varies rapidly near the edge of 
the gap, whose frequency we will denote by cjf,, which 
implies that the correlation time of the vacuum fluctua- 
tions can be comparable to the characteristic time scales 
of atoms, quantum dots or NV centers embedded in the 
PC structure. Hence, an accurate description of the dy- 
namics of systems in contact with highly structured envi- 
ronments, like the radiation field within PCs, may force 
us to go beyond the usual Markovian, and weak coupling 
approximations [TT|-[l3j. Therefore, a reliable method to 
describe the full system dynamics is highly desirable. 

Here, we demonstrate that the efficient, numerically 
accurate method developed in jTi - [l6j (TEDOPA: Time 
Evolving Density with Orthogonal Polynomials Algo- 
rithm) is an excellent candidate to describe the dynam- 
ics of quantum systems, such as impurity atoms and 



quantum dots, within structured reservoirs and quantum 
fields. To this end, we apply TEDOPA to describe the 
dynamics of a two level atom embedded in a photonic 
crystal. Our reason to chose such a system as a test bed 
is that, under certain limits and within the rotating wave 
approximation (RWA) , it is known to be exactly solvable, 
thus providing an excellent arena to prove the reliability 
of our technique. Most importantly, the method deployed 
here is valuable in describing the dynamics without the 
usual weak coupling and RWA. By evolving the system 
without invoking RWA, strong departures from the RWA 
result are found. We will show that when the system- 
environment coupling is strong, the atomic population 
dynamics differ from the RWA solution and an enhanced 
population trapping can be observed. When the system 
energy splitting. A, is much smaller than any other sys- 
tem and environment time scale, the RWA leads to an 
unphysical prediction for the system's coherence in the 
low frequency domain while TEDOPA allows to evalu- 
ate the accurate dynamics across the whole frequency 
domain. 



II. THE MODEL 

We consider a single atom or quantum dot, modeled as 
a two- level system with transition frequency A, coupled 
to the modified radiation field that exists within a pho- 
tonic crystal. The general interaction Hamiltonian can 
then be written in interaction picture as follows {h — 1) 
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where tr^ are the atomic spin ladder operators, /i = 1, 2 
are the two polarization modes of the light, and k is the 
field wave vector. In writing the Hamiltonian in this 
way, we ignore, for simplicity, any effects related to light 
polarisation and the existence of multiple bands in pho- 
tonic cystrals 0, H, [13 ■ The frequencies are defined as 
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ujk ± A, and we consider utk = Wf, + k'^/2m,, where of the total system can be written as, 



m is an effective mass acquired by a photon in a PC 
0, Q , corresponding to the dispersion relation of the ra- 
diation field for frequencies in the vicinity of a single 
band-gap edg e with frequency ujb, and k = |k|. Follow- 
ing Refs. |l7l f20l|. the microscopic form for the coupling 
between the light field and optical transition is given by 

fl'k.M = ^2fc,^e-^o'='/^ with nk,f, = A^J^^^^di2d ■ e^^^, 

where eg is the electric permittivity in vacuum, v is the 
quantization volume, d and di2 = XQke/2 are respec- 
tively the direction and magnitude of the dipole mo- 
ment with e the electric charge and Xq the Bohr ra- 
dius, and ek,;j for ^ — 1,2 represents the unitary po- 
larization vector corresponding to the two transverse po- 
larization modes. The factor e^'^ ^qI"^ physically arises 
from the fact that when calculating the coupling (see 
for instance [13] )j one should take into account that 
the electronic wave functions have a finite width Xq. 
Here, we have considered that the electronic wave func- 
tion for the electronic level j is approximately given by 



(r rj) /(2Xq)_ other words, the fac- 



tor e^^ ^a/'^ appearing in (7fc,p, often neglected within 
the dipolar approximation, is due to the fact that the 
electronic wave function varies on a scale that is approx- 
imately determined by the Bohr radius, Xq Consid- 
ering that the electron moves in a finite space leads to 
an effective momentum cut-off for the emitted photons, 
fco = 211 /Xq, and thus introduces a maximum photonic 
frequency ujq — uikg into the problem. 

The bosonic creation and annihilation operators of the 
photonic crystal light modes are 6k, and ^, respec- 
tively for a mode of wavenumber k and polarization fi. 
The above Hamiltonian is not exactly solvable and the 
RWA is usually considered, so that the fast rotating terms 
(i.e. terms ^ b\^^^a~ or ^j^^cr"*"), are neglected with re- 
spect to the other (energy conserving) terms. In this 

case, ff™ = Ek.M5k,p(^k,M'^"^'^''* + '*-^-)- ^hile 
the RWA is often a very good approximation, due to the 
existence of well separated time-scales, this is not always 
the case. In the case of strong coupling, we shall show 
that large deviations from RWA physics emerge, and that 
their description requires sophisticated numerical tech- 
niques. 



III. SINGLE ATOM DYNAMICS WITH THE 
RWA 

A single atom constitutes the simplest possible case, 
but represents a valuable example to prove the valid- 
ity or our method. As noted above, when consider- 
ing the RWA, the dynamics specified by the Hamilto- 
nian Eq.([T]) always remains in the one excitation sector 
and the time-dependent Schrodinger equation for the to- 
tal atom-photonic crystal system wave function can be 
solved exactly [1, [l^ • In this regard, the wave function 



\^{t)) = A{t)\i, {Ok}) -I- J2 BkAm Ik. 



(2) 



k,/j 



where |1, {Ok}) describes the atom excited and no photon 
in the radiation field, and |0, Ik,/^) represents no excita- 
tion in the atom, and a single photon in the mode k, /i. 
Note that this form for the system wave function is valid 
only when the RWA is assumed, and the dynamics is 
restricted to the one excitation sector. 

The time-dependent Schrodinger equation projected 
on the one-atom sector of the Hilbert space takes the 
form. 



dA{t) 
dt 

dt 
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Inserting the formal solution of the second equation into 
the first one, we have 



dA{t) 
dt 
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is the correlation function of the system-environment 
coupling [ll, [l^. We note that in order to calculate 
the correlation function, we have assumed that the fac- 
tor |d • ek.CTp/wfc changes very smoothly with k in the 
IBZ (first Brillouin zone) of the photonic crystal, and 
can be considered a constant factor in comparison with 



the rapidly oscillating exponential factor [17[. As dis- 
cussed in [l7j the validity of this approximation can 
be easily verified numerically, and it allows us to settle 
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Here, we have considered 



in addition that the sum over the two polarization modes 
gives rise to a factor two. Note that the high-frequency 
cut-off arising from the finite size of the electronic wave 
functions is required to avoid a singularity at the origin of 
time in the correlation function (such singularity appears 
for instance in j2l|). 

A change of variable between k and Uk, leads 
to expressing the correlation function as G{t) — 
— J(a;)e~*('^~'^^*, with the spectral density given 

by 
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where a = 27rr2^/(wQ^^). The amphtude A(t) can be ob- 
tained by taking the Laplace transform of Eq. ([5]) and 
inverting the resulting expression for the Laplace trans- 
form A{s) of A{t). From Eq. ([5]) we thus need to invert 
A{s) = A{0) / {s + G {s)) , where Q{s) is the Laplace trans- 
form of G{t). Once determined, the population of the 
atom in the excited state at time t is given by 

Using the Laplace transform method, a numerical 
solution for A{t) can always be obtained. In con- 
trast, an analytical solution requires additional approx- 
imations, namely that the system frequencies are very 
small compared to the cut-off frequency ujq, so that 
ujq ^ r^. A, a;f), and that all relevant time scales obey 
t :$> l/ujQ 0, [22l- Under these conditions, this leads to 
Goo{t) = _ae^('^-'^'')*+'^/'' 7^3/2, which is singular at the 
origin, but describes correctly times t >• 1 /wq 0, [H, [13 ■ 
Inverting the Laplace transform of A(^s) with this ap- 
proximation for G(t), the following solution is obtained 

m 

with /(a,AL,<) = (a/7r) \^ dx , ^"T-tr'^'l ■ Dcfin- 
ing r± = -(a/2) ± ^(a/2)2 - Al and Al = A - Wb, 
we have three different regimes: (i) If a^/2 > Al > 0, 
then ci = 0, (ii) If A^ > a^/2, then r\ — r_ and 
"^1 ^ 77^' ^'^^ finally, (iii) if Ai, < 0, ri = r+ and 
= Under these conditions, several dynamical 

regimes exist depending on the value of the parameter 
Al — Al — Ws, where is a frequency that arises from 
renormalisation of the optical transition by the environ- 
ment, ujs — 4r22/ti;o- Like the Lamb-shift encountered 
in traitional master equation approaches, Ws depends on 
the coupling strength and is discussed in Refs. [3 I3| 
and in Section |Vl In this regard, for Al > the atom 
is not excited in the stationary state (i.e. relaxes com- 
pletely to the ground state), whereas for negative values 
of Al there is a stationary residual excited atomic pop- 
ulation (corresponding to a situation where the emitted 
photon remains localized nearby the atom). Thus, tuning 
A L from negative to positive values leads to a cross-over 
between two distinct regimes [13, [l^ . 

A. Chain mapping, rotating wave dynamics and 
TEDOPA 

When the RWA is not applicable, the dynamics of 
the global atom-light system cannot be described in the 
one-excitation manifold, and the direct wave function 
approach presented above becomes intractable. How- 
ever, a new technique has recently emerged in the the- 
ory of open quantum systems which allow this task to 
be performed with numerical exactitude. This tech- 
nique [3, [l5| employs time-adaptive density matrix re- 
normalisation group (t-DMRG) algorithm to compute 



the evolution of the full atom and light field wavefunc- 
tion. The t-DMRG algorithm enables extremely accu- 
rate simulation of large many-body systems with nearest- 
neighbour interactions in 1_D, and has become widely 
used in cold atom and condensed matter theory (25j . 

To apply this technique to the present problem, we 
start with the full atom-light field Hamiltonian in the 
continuum Schrodinger picture, 

H = ^{l + (7,)+a^ f dkh{k){h{k)+h{k)'^) 
2 Jo 

+ [ uj{k)b{k)^b{k) (9) 
Jo 

where the field operators obey that [b{k),b{k'y] = S{k — 
k'), h{k) is the atom-light coupling strength and cj(fc) is 
the dispersion of the photons. The parameter k labels the 
modes and is normalised so that the highest frequency of 
the modes, that we will denote by uic, occurs at fc = 1. 
When the environment is initially in a gaussian state, it 
is possible to integrate out the environment degrees of 
freedom and derive a formally-exact expression for the 
reduced state dynamics of the atomic system using path 
integral methods [ssj . Importantly for the application of 
t-DMRG, this result shows that the effects of the pho- 
tonic environment will only enter through the spectral 
function J(w) of the photons [33]. In the following sim- 
ulations, we assume that the initial state of the total 
atom and environment is a product state of an atomic 
state and the zero-temperature vacuum state of the light 
field. As we are only interested in the reduced state of 
the atom, we can make use of the fact that the only rel- 
evant environmental property is the spectral function to 
represent the bosonic environment in a form that allows 
us to transform it into an equivalent system suitable for 
TEDOPA simulation. 

In terms of the continuum dispersion and coupling 
functions, the spectral function defined in Eq. ^ takes 

the general form J{uj) — h'^{t{uj)) '^'^^'^^ , where e{k) is the 
inverse function of the dispersion uj{k), i.e. e{uj{x)) = x 
[Til l . From this one can see that a variety of different pairs 
of h{k) and uj{k) can lead to the same spectral function, 
and we exploit this. Taking h{k)'^ — LdcJ(<^{k)) and the 
linear dispersion Lu{k) = uib + ujck, with ojc the max- 
imum frequency of the environment considered in the 
simulations, it can be found that the spectral function 
characteristic of ([9]) (as described in [1^) is given by 
J{uj)Q{uJc — where J{uj) is the spectral function of 
Eq.(l7l) and Q{x) is the Heaviside function. As the dy- 
namics of the reduced atomic state only depend on J{uj), 
the t-DMRG simulations (with these choices of aj(A:) and 
h(ky) simulate exactly the same reduced dynamics im- 
plied by Eq. ([ij with the physcial microscopic dispersions 
and coupling functions described in Section |TT1 Because 
of the hard cut-off in the spectral function, dynamics 
are captured accurately on all frequency scales up to Wc, 
which is chosen to be much larger than the physcial scale 
loq that determines the frequency range over which J{lij) 
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is significantly different from zero. Therefore, by tak- 
ing Wc ^ the t-DMRG describes ah relevant physics 
beyond times of ss 1 /w^ with numerical exactitude. 

We then perform a unitary transformation of the 
Hamiltonian given in Eq. ^ by defining new en- 
vironment modes a„ = p^^ h{k)'Kn{k)b{k), where 
7r„(fc) are monic orthogonal polynomials which obey 

dk J{g{k))TTn{k)-Km{k) = PnSnm- By choosiug a liucar- 
in-k dispersion, the properties of the monic orthogonal 
polynomials ensure that this transformation is both real 
orthogonal and leads to a nearest-neighbour chain Hamil- 
tonian given by [l^ [13 , 

Hchain ^ ^i'i^ + o-z) + gcTxiao + al) (10) 

oo 

Within this picture, the two level atom now couples to 
one end of a chain of coupled harmonic oscillators. The 
frequencies e„ and couplings i„ are simply determined by 
the recurrence coefficients of the monic orthogonal poly- 
nomials which are effectively determined by the spec- 
tral function (in the language of orthogonal polynomi- 
als, the spectral function is treated as a weight function) 
[ISL 1231 . The use of the theory of orthogonal polynomials 
also allows rigourous proof of this mapping for almost 
any physical spectral function, including continuous, dis- 
crete, mixed and gapped spectra [H, H^] . We also note 
more generally that this procedure effectively enables a 
highly efficient way of simulating one-dimensional quan- 
tum field dynamics coupled to non-linear objects, such 
as spins, and also other quantum fields. 

After this procedure it now becomes possible to ap- 
ply time-adaptive density matrix renormalisation meth- 
ods J3) HE] to simulate the unitary evolution of the to- 
tal atom and chain wavefunction |5'(i)). From |^(i)), 
the real-time expectation of any observable Oa of the 
atom can be simply calculated by evaluating {Oa{t)) = 
{'i'{t)\Oa\'^{t)). Although we will not pursue this here, 
another advantage of this method is that we also have 
full access to the evolving state of the environment field 
itself. This can provide deep insight into real-time dissi- 
pative processes, as has been recently demonstrated in a 
study of exciton transport in photosynthetic complexes 
[2^ . The combination of the powerful analtyical mapping 
with the numerical strength of t-DMRG we shall hence- 
forth refer to as the TEDOPA (Time Evolving Density 
with Orthogonal Polynomials Algorithm) technique, and 
a comprehensive overview of its theory, implementation 
and applications can be found in p7| . 

B. A quick comparison of methods 

While numerical simulations with t-DMRG can be 
made arbitrarily precise in principle, in practice fi- 
nite bond-dimension employed in the simulation will 



limit precision. The numerical simulation error can be 
bounded from above but the resulting bounds are very 
conservative as they are growing in time. Hence it is 
instructive to compare simulation results against know 
exact solutions. 

To this end, we consider the RWA approximation 
and compare the results obtained using TEDOPA and 
the numerical solution of A{t). As the analytical (in- 
verse Laplace transform) solution for A{t) discussed in 
the previous section is only valid when t ^ 1/wo and 
luq D,,A,LUb, we use the Piessens method [23| to nu- 
merically invert the Laplace transform of A{t) for the 
large but finite values of ujq used in our TEDOPA cal- 
culations. In addition, it is possible to numerically inte- 
grate the RWA Heisenberg equations of motion in the 
chain representation of Eq. (|lip to obtain a numeri- 
cally accurate evolution of the atomic population. Fig- 
ure [T] shows the population (|A(t)p) in the excited state 
of the atom for a = 1, wq = 100, gap-edge frequency 
cjf, = 5, cjc = 800 and A = 1 calculated by t-DMRG, 
the Piessens method and the numerical integration of 
the Heisenberg equations. The figure shows that the 
t-DMRG result converges to both exact RWA results, 
while the inset shows that the finite system-size recur- 
rence features which appear in the Heisenberg integration 
scheme are reproduced almost exactly if the same num- 
ber of sites are used in the t-DMRG simulation. We note 
that it is also possible to perform t-DMRG in the Heisen- 
berg picture, and that for the RWA example considered 
here it will be numerically accuarate [29 . 30] . Anticipat- 
ing later results, which show that TEDOPA captures the 
non-rotating wave physics extremely well, we might also 
expect Heisenberg t-DMRG to perform well in the non- 
rotating wave case [2^, [30| . The slight departure between 
the t-DMRG and Heisenberg method from the Piessens 
method is due to the growing inaccuracy and instabil- 
ity of the Piessens method (a general feature of many 
numerical Laplace inversion techniques) as time evolves 



IV. ACCURATE NUMERICAL RESULTS 

The cross-over occurring when tuning A/, from nega- 
tive to positive values can also be observed when con- 
sidering the exact numerical solution of A(t), but until 
now it had only been described within the RWA. Thus, 
a question that remained open is whether this transition 
is a consequence of the RWA, or can also be observed 
in situations where the RWA is invalid. To explore this 
we now present results in the regime of strong coupling 
where we expect, and shall demonstrate, that deviations 
between RWA and the full simulations to be most appar- 
ent. This regime of strong coupling is reached when the 
energy shift induced by the environment coupling when 
A = is comparable to the energy splitting of the optical 
transition when A is finite. The environmental energy 
(polaron) shift is Een = J{^)^~^ ~ ce{^o/'^)'^ 
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FIG. 1. (Color online) Comparison between t-DMRG (solid 
black), Heisenberg (dashed red), and Piessens (blue dots) so- 
lutions for |j4(t)p, for the parameters a — 1, u>o = 100, gap- 
edge frequency ui, = 5, cUc = 800 and A = 1. Time is in units 
of 




FIG. 3. (Color online) Stationary state atomic population for 
the RWA (red squares) and the complete Hamiltonian (blue 
dots) as a function of A. The initial state is a fully excited 
system, and we have used a = 1, ujo — 100, ujt = 5, and 
= 800. 
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FIG. 2. (Color online) Excited state population dynamics for 
the RWA (filled plot markers) and full Hamiltonian (empty 
plot markers). For these data a — l,iuo — 100, ujb = 5, and 
LOc = 800. Time is in units of 1/a^ 



|16j . and by taking parameters a = 1, wq = 100, cj;, = 5 
in the rest of this article, we ensure that we remain be- 
tween the strong and intermediate couphng regimes for 
aU of the values of A we shall consider. 

Fig. [2] shows the time evolution of the population in 
the excited state as a function of time, when starting 
from the excited atomic state, and for fixed bath couphng 
and variable A. Clear differences in the evolution can be 
seen, particularly in the final values of the excited state 
population. Indeed, when A = 15 the result of the RWA 
predicts no population in the excited state, whereas the 
dynamics for the full Hamiltonian ([1]) still shows a non- 
zero final excited state population. In general, the RWA 
overestimates the excited state stationary population for 
Ai < 0, and underestimates this quantity for A^ > 0, as 
seen for A = 10. We also note the oscillatory features of 
the population dynamics, which in both cases decrease 
in frequency as A approaches the band edge. However, 



at fixed A the RWA and full Hamiltonian frequencies 
are markedly different, especially when A — ujf,. This 
behaviour will be rationalised in section [V] 

Fig. [3] shows the stationary atomic population for an 
initially excited atom, with respect to the atomic fre- 
quency A. In this Figure, the mentioned cross-over from 
a photon-atom bound state to a total relaxed state is 
sharper within the RWA than when the full Hamiltonian 
Eq. ([T]) is simulated. A clear difference in the residual 
populations above the cross over is also observed. As 
we shall see, these differences arise due to the improper 
treatment of the renormalisation of the atomic transi- 
tion energy in the RWA. However, at very weak cou- 
pling, these differences are negligible and the cross-over 
happens at A « Wf, (not shown). We also note that a 
sharp 'transition' at A^ =0 only appears as — > oo 
[1, [l3|, and when the coupling is very weak. However 
one should also note that the cross-over width in Fig. [3] 
is over-emphasised by the logarithmic scale. 

FiglS] also shows that the residual population in the 
excited state remains larger in the simulations of the full 
Hamiltonian relative to results obtained within the RWA 
approximation. This arises due to multi-photon effects 
which cannot be described in RWA. For A ^ w^, the 
atom is able to exchange energy with the light field via 
photon emission/absorption, and at long times relaxes 
towards the ground state of the full atom-photon Hamil- 
tonian. However, the interactions with the bath alter the 
effective eigenstates of the system, leading to the residual 
excited state populations in the atom-light field ground 
state. This state can be approximately described by the 
variational polaron ground state, first suggested in [sif . 
Following the procedure set out in [3l| , the reduced den- 
sity matrix of the atom's variational polaron ground state 
- in the basis |±), where |±) are the eigenstates of ax - is 
givenbyp^*°" = ^l+>(+l + ^|->(-|-^|-)(+|-*l+)(-| 
where $ = A/A and A is determined from the implicit 
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equation 

r 2 /"^ 
A = A exp / 



The ofF-diagonal elements of the density matrix are 
suppressed by dressing, or polaronic, correlations (state- 
dependent coherent displacements) between the bath and 
atomic states |±). These correlations effectively suppress 
the effective energy gap (A — > 0) between the atomic 
states I t) and | J,) through the reduced overlap of the dis- 
placed light mode wave functions which dress the states 
|±) [ll,|3l|. From the density matrix yOg*"™, the popu- 
lation in the excited atomic level is = 1(1 - For 
large a, w;, ^ A <C the approximate self-consistent 
solution of Eq. ^ is A w This shows 

that as A gets larger, the re-normalisation from the en- 
vironment gets smaller, so $ — ?> 1 and P-|- — >■ 0. For the 
parameters in Fig. [3l we find that the solution of Eq. 
dTTI) with A = 30 gives = 0.026 which agrees with the 
data very well. Note that these formulas are only valid 
when the atom is able to relax to the polaronic ground 
state, which occurs for A > Wf, [13]; at A < Wb, relaxation 
is blocked by energy conservation. In this case a better 
prediction for the residual population is to assume that 
the final state is the dressed excited state, in which case 
the final population is = ^(1 -I- with $ obtained 
from Eq. (|lip . In the RWA approximation the analytical 
theory predicts a total de-excitation of the atom at large 
A > Wb. This is seen in Fig. [3l although the transition 
is slightly smoothed out by the finite value of wq used in 
the simulations, as discussed above. 
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FIG. 4. (Color online) Time evolution of the absolute values 
of the coherence {(Jxit)) for RWA (filled plot markers) and 
the full Hamiltonian (empty plot markers). For these data 
a — l,cJo ~ too, iUb = 5, and ujc — 800. Time is in units of 
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FIG. 5. (Color online) Frequency of coherence oscillations 
using the RWA (red squares) and the complete Hamiltonian 
(blue dots) as a function of A for an initially-prepared super- 
position of ground and excited states. Parameters are a = 1, 
Ldo ~ 100, ujb ~ 5, and uuc = 800. 



V. COHERENCE DYNAMICS AND THE 
FAILURE OF THE RWA 

Let us now consider the evolution of the system co- 
herences {<Tx{t)), starting from an initial condition where 
(o'2;(0)) = 1. In the absence of the radiation field, we 
would expect that {crx{t)} — cos(At) and in its presence 
these coherent oscillations would - in a simple markovian 
picture - become damped and shifted to lower frequen- 
cies [nl, nil- Figure m shows the dynamics of \{ax{t))\. 
The difference in dynamical behaviour is dramatic. As 
A is increased from well below the gap to well above, the 
full simulation shows monotonically increasing frequency 
for the coherence oscillations, with some very weak non- 
markovian dephasing just below the band edge and much 
faster, exponential decay for A = 15. In stark contrast, 
the RWA oscillations show a non-monotonic dependence 
on A, initally decreasing in frequency and then increasing 
as A increases. 

Figure [S] shows the frequency of the coherence oscilla- 
tions for the full Hamiltonian and the RWA for different 
values of the atomic frequency A. The behaviour seen for 
the RWA curve in Figs. H] and [S] can be understood by 
looking at the pole structure of the exact RWA solution 
expressed as a Laplace transform. For our initial condi- 
tion, the Laplace transform of {ax{t)) can be expressed 
as 

'C[(fT,(t))] = i(s + iA + C[G{t)]y' + h.c. 

The poles oi C[{ax{t))] determine the dynamics of the co- 
herence, with their imaginary parts giving the oscillation 
frequency and their real parts quantifying the dephas- 
ing rate. Although the Laplace transform C[G{t)\ can be 
computed exactly, we can qualitatively understand the 
key features of Fig. [5] by considering just an expansion 
to lowest order in the small quantities A/ojo, ^b/'^a- This 
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gives 



C[G{t)]) — iaW h aV is — cjf,, 

TT 



with a given in Eq. When A 

C[{ax{t))], s±, are approximately 



0, the poles of 



s± « zLia{i 



This shows that the oscillation frequency converges to 
a finite value as A — 0, as is observed in Fig. [Sj The 
poles are completely imaginary, and these oscillations are 
undamped, although there is a slight loss of amplitude 
due to dynamical transients. For small but finite A, the 
poles now appear at approximately 

s± « ±i(A — a\ — + Q!V^A~~Wb), 

V TT 

leading to a reduction of the effective oscillation fre- 
quency as the A term compensates the frequency shift 
coming from the environment. The large dip in oscilla- 
tion frequency seen in Fig. [3] corresponds to a significant 
cancellation between these terms. For large values of A 
(A ^ ujh), the poles are dominated by the contribution 
from A and are approximately 

s± w ±iA + 

The oscillation frequency is now proportional to A and 
has acquired a real part which damps the oscillations with 
a rate given by the Fermi Golden Rule as F = J(A). This 
is the result one would obtain by a standard weak cou- 
pling and Markovian approximation to the open-system 
dynamics with the spectral density of Eq. ([7]). 

The behaviour of the oscillation frequency for small 
A is clearly a pathology of the RWA approximation and 
arises because of a poor treatment of the low frequency 
modes in the problem. The breakdown of the RWA man- 
ifests itself in the renormalisation of the atomic energy 
A which leads to a shift in A that is larger than A it- 
self, an effect which produces the minimum in Fig. [5] 
and the convergence to a finite oscillation frequency as 
A — > 0. The correct behaviour in this limit can be calcu- 
lated exactly for the full Hamiltonian ([IJ, as for A = 
the model corresponds to the exactly solvable indepen- 
dent boson model [s^. This shows, as one would expect, 
that {ax{t)) does not oscillate and remains constant at 
A = 0. This limiting behaviour is correctly described 
by the TEDOPA simulation of the full Hamiltonian. For 
small, but finite A, the frequency simply increases mono- 
tonically with A and the TEDOPA results agree very well 
with the prediction 

A« Aexp(-a/V^) + 0(^^) 

given by adiabatic renormalisation theory (a simplified 
version of Silbey-Harris theory which is valid in the small 
A limit) fQ. 



We finally remark that the failure of the RWA co- 
herences also manifests itself in the frequencies of the 
populations we presented in Fig. [21 Neglecting dissipa- 
tive effects, the oscillation frequency in the excited state 
population below the band gap would be expected to 
vary as Al — (A — Wb), and would thus vanish as the 
renormalised atomic energy moves into the bandwidth 
of the light field. By comparison with Fig. [SJ it can be 
seen that for the values of A considered in Fig. [2l the 
over-suppression of A for the RWA leads to both a sup- 
pression of population decay (atomic transition is further 
away from the band) and, consequently, spuriously fast 
population oscillations. This is best seen in Fig. [2] at 
A = cjb = 5, where the oscillatory part of the RWA pop- 
ulation dynamics is almost four times faster than in the 
full Hamiltonian simulation and retains 50% more of the 
excited state population. This failure of the RWA renor- 
malisation for low values of A also leads to artificially 
large residual population in the excited (relative to the 
full simulation) seen in Fig. [3] for A =< 11. 



VI. APPLICATION TO PHOTONIC CRYSTALS 
AND STRONGLY-INTERACTING LIGHT 
MATTER SYSTEMS 

Considering a photonic crystal with a gap in the opti- 
cal region, couplings of the order of = 10^ — lO^Hz 
may be achieved, which in our units [35] would be 
= 10~^ — 10~^°. In order to enlarge the cou- 
pling, and thus reach the regime discussed in this pa- 
per, there are different possibilities. One may consider 
for instance photonic crystals with narrow band-widths 
Aoj, what would give rise to smaller cut-off frequencies 
Wo = uJko = + ZAwfeg/S. Here we have considered the 
fact that the h/2m = PAoj/3 [i3|. 

Strong system-environment coupling can also be 
achieved when strong collective effects produce an en- 
hanced effective system-environment coupling [sgI ISTj . 
The simplest situation in which this enhancement can 
be achieved is by considering a collection of N atoms 
uniformly interacting with the radiation field or environ- 
ment, or similarly, when describing the collective cou- 
pling of a spin system to superconducting resonators [s^ . 
In this case, for a single excitation in the system (i.e. 
only one atom in the ensemble is excited), the situation 
is completely analogous to the one described here, ex- 
cept for the fact that the coupling is enhanced by a fac- 
tor in the photonic crystal case 0- Thus, a value 



^off 



1 would be achieved if we have a high 



enough number of impurity atoms. Of particular inter- 
est is also the multi-excitation regime for this configura- 
tion. A spontaneous polarization phenomenon, in which 
the atomic coherences grows from zero to a finite value, is 
predicted to occur for atoms in photonic crystal-like envi- 
ronments, when the atomic frequencies are resonant with 
the band edge A = and when considering the semi- 
classical and RW approximations 0, [IB [3 ■ The results 
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in our work show strong differences in the behaviour of 
single atom coherences with and without considering the 
RWA, particularly for small Ar , and suggest that apply- 
ing the novel technique [13, ^M, for treating the system- 
environment interaction in the analysis of multi-atomic 
systems will shed new light into this and other intriguing 
collective phenomena. Other possible enhancement fac- 
tors for the coupling coefficient would involve placing the 
atom within a photonic crystal cavity [s^ . In a more gen- 
eral scenario, one may consider a situation where the ef- 
fective system-environment interaction, characterised by 
y/rjQ in our chain representation of the system, is larger 
than the typical frequencies of the light field, ujb- As 
we have shown, in this regime the ground and excited 
atomic states are entangled with the light field through 
polaronic correlations, and these are very similar to the 
shifted vacuum states which appear in the recent theory 
of ultrastrong coupling in quantum wells and circuit QED 
where this regime can be realised [40, 41]. In conclusion. 



this work has demonstrated that our TEDOPA method 
for the treatment of system environment interaction has 
allowed us to reveal strong deviations between the exact 
and approximate solutions in physical models of practical 
importance, and hence suggests its use in a wide variety 
of settings where non-perturbative system environment 
interactions are expected to play a role. 
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